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Abstract 

In this paper, based on a discrete total variation model, a modified dis- 
cretization of total variation (TV) is introduced for image processing prob- 
lems. Two optimization problems corresponding to compressed sensing mag- 
netic resonance imaging (MRI) data reconstruction problem and image de- 
noising are proposed. In the proposed method, instead of applying isotropic 
TV whose gradient field is a two directions vector, a four directions dis- 
cretization with some modification is applied for the inverse problems. A 
dual formulation for the proposed T’V is explained and an efficient primal- 
dual algorithm is employed to solve the problem. Some important image 
test problems in MRI and image denoising problems are considered in the 
numerical experiments. We compare our model with the state of the art 
methods. 
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1 Introduction 


There are different methods to solve image processing problems such as spa- 
cial filtering {[11, 6], methods based on partial differential equations [21, 26], 


*Corresponding author 
Received 1 January 2020; revised 1 March 2020; accepted 31 May 2020 


Alireza Hosseini 
School of Mathematics, Statistics and Computer Science, College of Science, University of 
Tehran, Iran. e-mail: hosseini.alirezaQut.ac.ir 


Erfan Ebrahim Esfahani 
School of Mathematics, Statistics and Computer Science, College of Science, University of 
Tehran, Iran. e-mail: ebrahim.esfahani@ut.ac.ir 


+ This article was suggested by the Scientific Committee of the Third Iranian Sem- 


inar on Control and Optimization for publication in IJNAO, which was accepted after 
independent review. 


87 


88 A.R. Hosseini and E.E. Esfahani 


variational methods [8, 2, 10, 1, 5], and machine learning approach as well as 
the methods based on deep learning [17, 25]. Variational methods propose 
some different functionals as the tools for calculating accidental variations in 
a signal or function. The pioneer work in this area was proposed by Rudin, 
Osher, and Fatemi [24]. Recently, Condat [10] proposed an isotropic discrete 
total variation with high accuracy. Moreover, the total generalized variation 
as a generalization of total variation functional has been proposed in [5]. 
Consider the following optimization problem to solve a continuous version of 
mathematical image problems: 


. 1 2 
min [ |G(s(t)) — g(t)Pat + J(s), 


where G is an appropriate function (for more details see [7]) and 
J(s) = sup i-f s.divodt : ¢ € C1(Q, RY), |o(t)| < 1, Vt € a} i 
Q 


Also J is the dual definition of total variation (TV) of the £1(Q) (locally 
Lebesgue integrable) function s. If G = I, then problem (1) is named denois- 
ing problem (for image g). A function s is said to have bounded variation 
whenever J(s) < oo. The space BV(Q) of functions with bounded variation 
is the set of functions s € L1(Q) such that J(s) < oo, endowed with the norm 
IIslleva) = IIsllz2@@) + J(s). Obviously, for the smooth function s € C1(Q) 
(or se W1(Q)), 


ioe i. \Vs|dt 


the minimization of J(s) is equivalent to minimization of the sum of abrupt 
variations (the majority of the derivative) over its dimension. To solve real 
digital imaging problems, discretization of (1) is inevitable. There are differ- 
ent forms of discretization for TV functional in the papers. As the review of 
such literatures is long some, we refer the readers to [10, 13] and references 
therein. 

Total generalized variation (TGV) is a cleverly defined functional for cal- 
culating a combinational value of function variation by means of derivatives 
up to a given order. TGV has some superiorities in comparison with the 
first-order TV, such as its ability to distinguish staircase artifacts, and as 
a result, more cleaner image can be reconstructed. Usually in applications, 
the second-order TGV is sufficient, because the higher-order will face to the 
complexity of calculation. The second-order TGV is defined by the following 
optimization problem (see [5]): 


TGV2(u) = sup iff u.div?udx : v € C2(Q, Sym?(R*)), ||div' vlloo < az, = 0, i} : 
Q 
In [5], a discretization of the TGV was proposed. Based on discretization of 


TV in [13], in this paper, a new discretization of TV (1) is given; operators 
are changed and linear operator QO, is inserted to ensure the projection on the 
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vertices. Then the magnetic resonance imaging (MRI) data reconstruction 
and image denoising problems will be formulated based on the new proposed 
discrete TV and will be compared with some state of the art models. In the 
proposed model, instead of using two directions gradient field in isotropic TV, 
a four directions discretization model is proposed. Actually, we can take more 
information from structure of the images, and consequently a new function 
for measuring the majority of the noise and oscillation is proposed. Moreover, 
after some modification based on the Condat’s model [10], a discretization of 
the total variation with better isotropy property will be defined. Numerical 
results show that the proposed method has better performance in comparison 
with isotropic TV, TGV, and some other famous models in terms of removing 
noise and compressed sensing MRI problems. 

The paper is arranged as the following; in Section 2, the Chambolle—Pock 
algorithm is explained for solving optimization problems in image processing. 
In Section 3, a new discrete TV model is proposed. In Section 4, the problem 
of reconstruction of undersampled MRI data as well as the problem of image 
denoising are formulated. Some important test problems in medical MRI and 
denoising are solved numerically and compared with some state of the art, in 
Section 4. Finally some conclusions are given in Section 5. 


2 A numerical algorithm for solving convex image 
problems 


2.1 Algorithm description 


Consider a finite-dimensional vector space X with an inner product (-,-)x 
and Euclidean norm ||-||2 = \/(-,:)x. C C X we define the indicator function 


of C as 
0 rEC 
6 — 7 ? 
o(@) e cE C. 


For any function F' : X — [—oo,+00], the convex conjugate at z € X is 
defined by 
F*(z) = max (a, z) — F(a). (2) 


Moreover, for any 0 > 0 and & € X, the proximal (or proximity) mapping of 
F at Z is defined by 


Ile = =|)? 


30 + F(z). 


proxy -(Z) = arg min 
x 


Now, suppose that Y is another real vector space with inner product (-,-)y 
and the corresponding induced norm |] - ||2 = \/(-,:)y. Let K: X — Y bea 


90 A.R. Hosseini and E.E. Esfahani 


bounded linear operator with the operator norm 
|||] := max{||Ka||2: a2 € X and ||2||2 < 1} < co. 


The adjoint of K is defined as an operator K* for which the equality 
(Ka,y)y = (a, K*y)x holds for all x and y. The general problem, we tackle, 


is 
min F( Kx) + G(x), (3) 


where F' : X — [0,00) and G: Y + (0,00) are proper, closed, and convex. 
By virtue of (2), the convex minimization problem (3) can be reformulated 
as the saddle point convex-concave problem 


min max (K2,y)y + G(x) — F*(y), (4) 

coy 
which is also known as the primal-dual formulation of (3). Hereafter we 
assume that a solution to this problem always exists. The method we employ 
to solve (4), which lays the groundwork of what we shall propose in this paper, 
is the celebrated Chambolle—Pock algorithm, recalled below: The algorithm 


Data: K, k*,F*,G 
Result: x* 
initialization: Choose parameters ¢ > 0, 7 > 0, and initial estimates 
(2°, y)EXxY. 
while convergence criterion not met, for k=0,1,..., do 

L-y*** = prox, p.(y* + 0K"), 

2. ¢°1! = prox_p(a* — 77K *ty**), 

Sehr! = Det tt a. 
end 

Algorithm 1: The Chambolle—Pock algorithm for solving problem (4) 


is guaranteed to converge provided that oT < TIP: In practical situations, 
where computing the exact value of || K|| is not easy, finding an upper bound 
B>||K||? and setting o =7 = Ta are sufficient for convergence. We chose 
this algorithm over other well-known first-order methods such as ADMM 
or FISTA because of its superior convergence rate [9], and the fact that it 
requires fewer parameters to be tuned. 


2.2 Convergence analysis 


Here, we express the convergence theorem of the Chambolle—Pock algorithm 
to solve problem (4). The convergence analysis was argued in [9]. The fol- 
lowing theorem guarantees the convergence of Chambolle—Pock algorithm for 
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solving problem (4). 


Theorem 1. [9] Assume X and Y are two finite-dimensional Hilbert spaces 
and that T,0 > 0. Moreover assume o7||K||? < 1. Then there exists a solution 
(Z,y) € X x Y of problem (4) such that the sequence (x”, y”) in Algorithm 
1 converges to (Z, 9). 


Note that, in [9], it was shown that the rate of convergence is O(+). More- 
over, under some uniformly convexity assumptions, the rate of convergence 
can be improved to O(=y). 


3 A Discretization for TV 


Assume that x € R%'*? is a gray scaled two-dimensional image and that 
A= {(n1,n2),n1 =1,...,Ni,n2 =1,..., No} is a set of grids. Consider the 
following new discrete TV model, which is based on four directional gradient 
with imposed some operational constraints: 


TVais(%) = Max, cgaymixre{< Dz,u >: |Ozu(ni,n2)| <1, [Onu(ni,n2)| < 1, 
|O.u(ni,n2)| <1, |Ozu(ni,n2)| < 1,V(m1,n2) € X}, 


where 
D= (D,, D2, D3, D4), 
Dyx(n1,n2) = x(n, + 1,n2) — x(n1, n2), 
Dox(n1,n2) = x(n, ng + 1) — x(n41, n2), 
D3x(n1,n2) = «(ny +1,n2 +1) — x(n1, 2), 
) = a(n — 1,n2 +1) — x(n1, no), 


4 
|\O,u(n4,2)| = S“[(Oxu); (ma, na)? (n1,n2) € X,* =f,0,¢,4, 


and 

(Oyu)1 (M1, n2) = ui(n1, ng), 

(Oyu)2(n1, n2) = + (u2(n1, n2) TT ug(n1, N2 = 1) 
+up(n1 + 1,2) + u2(n1 +1,n2—-1)), 
(Oxu)3 (m1, M2) = 5 (U3 (M1, N2) + U3(M1,N2 — 1)), 
(Osu)a(m1, 2) = Z(Ua(M1, M2) + U4a(M41, M2 — 1) 
bua (ny +1, n2) + ua(n41 +1,n2—- 1), 


wi 


Joon 
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ui(1, N2) + ui(n1 — 1, nz) 
1) t ur(n41 1, ne t 1)), 


+ Ale 
S 
an 
3 
es 
= 
1) 


(O.,u)2(n1, n2) = ua(n1, n2), 
(O.,u)3(N1, n2) = 4 u3(n ,n2) ar u3(n1 — 1,n2)), 
(O.,u)a(mi, n2) = $(ua(m1, n2) + a(n + 1,12), 


uz (M1, M2) + u1(N41, M2 + 1)), 
ug(n ,n2) ai U2(n1 + 1,n2)), 


aa ug(n1, N2 — 1)), 
(u3(m1, 22) + ug(n1,N2 — 1) 

u3(n1 1,n2) t u3(n41 1,n2 1)), 
et)a(m1, N2) = 4(u4(n1, n2) + ua(ny + 1,n2) 
+u4(n1,N2 = 1) + ua(n41 + 1,n2 aa 1)). 


= eae + ui(m1 —1,n2)), 
f 


Theorem 2. Assume that F' is convex and that \ > 0. Then the problem 


min F'(v) + ATVais (2), 


is equivalent to the following optimization problem: 


min F(a) + A{|vz| + |ves| + |ve| + joy}, 
subject to 
OF (vp) (m1, n2) + OF, (Yer )(m1, M2) + OF (ve) (M1, n2) 
+04 (v+)(m1,n2) = Dx(m1,n2), (mi, n2) € X, 


where 


v¢(n1, 2) 


g(¢ (M1, na) + UF (M1, M2 +1)+ Up (M1 —1,n) 
ug (na 1, ne 1) 


OF (vy) (m1, 22) = 5 
3 (U¢(n ,M2) + Us (ni, M2 +1)) 


(s(n, ng) + vs(na, mg+1)+ up(na —1,n2) 
up ( 1,n2 + 1)) 
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4 (v2, (1, m2) + ve, (m1 +1, 2) + v2 (m1, m2 — 1) 
+l, (ni + 1,2 —-1)) 


2 
O%, (ve) (1, N2) = uz, (M1, N2) . 


$(v3, (m1, N2) +v3,(n, + 1,n2)) 


$ (ve, (m1, 22) + vé,(n, — 1,7n2)) 


$(v2(n1,N2) + v2(n1,n2 + 1)) 


O% (v2) (m1, n2) = + (v3 (ni, n2) + u3(n1,n2 +1) + v3(n1 +1,n2) 

+u3(n1 + 1, neg +r 1)) 

7 (ve(m1,n2) + ve(n1 — 1,n2) + ve(n1, m2 + 1) 
+vue(n1 = 1, neg +r 1)) 


L(y (m,n2) + vi (m,n — 1)) 


$ (v4. (ni, n2) + 4(n1 — 1,Nn2)) 
O7 (v4) (m1, n2) — 7 , 
Vi(N1, 12 


vt (m1, n2) 


|v.| = S- [v2 |?,* =t, 0, ¢, +. 
n1,N2 


Proof. It is a direct conclusion of Fenchel—Rockefeller, primal-dual theorem 
(see [10], for instance). 


and 


4 Formulation of some image processing problems 


4.1 Compressing MRI data 


Consider the following problem for reconstructing MR images from under- 
sampled k-space data: 


. tL 
min 5<||Fmu — 13 + Bl|Wullr + TVais(u) (5) 


94 A.R. Hosseini and E.E. Esfahani 


where u € U = R”*” is the MR image, which is defined in spatial domain and 
Fm is the undersampled Fourier operator, which converts signals in spatial 
domain to its corresponding frequency representation and undersample it. In 
other words, it can be defined as F,,u := M © F(u), where F is the two- 
dimensional fast Fourier transform, Mx, is the k-space sampling mask with 
ones at m sampled frequencies and zeros at undersampled locations, and © 
is the component-wise multiplication. We further assume that the sampling 
rate ™; is considerably small (highly undersampled k-space) and that b € 
C”*” is the partially scanned k-space with only m sampled frequencies and 
n? —m zeros. The so-called zero-filled solution is defined as tgp := F—1(b). 
Furthermore, W is a discrete wavelet transform, and the constants 6, \ € R 
are regularization parameters. 

Referring to Theorem 2, we can get the following equivalent formulation: 


min 5h | Fme b|3 + Bl Mulla + lotls + osha + los + leds, 
s.t. OF(vy4)(m1, n2) + OF, (v.,)(n1, N2) + OF(ve) (M1, n2) (6) 


+ O% (v4) (1, n2) = Du(n1, 2), (n1, n2) EA. 


A= {(n1,n2):ni,ng =1,...,n}. 


To solve this problem with Algorithm 1, we need to reformulate problem (6) 
to the form of (3). Now we define the following linear operator K: 


u € RM XNe2 


vw 0 0 0 0 vy € (R*)™ x No ke (IR)N1* Ne 
Ke=|Fm0oO 0 0 0 Vex € (IR*) N12 he (IR)N1x No _ 
“DOL O% Of OF J | v6 (RO BE (RA) MxNe 


U4 E (Rare Ne 


Assume 1 
FY) = sy ll¥2 — bI5 + BULYilla + 540} (¥3), 


G(X) = Jogla + lud|a + lochs + loys. 


From arguments about calculating adjoint of a given operator, we get 


: r 
F*(y) = 5 llyell2+ < y2,b > +6).)..<e(yi)+ < y3,Du >. 


4.1.1 The y-subproblems 


According to Algorithm 1, for {Y} = (k,h,w)', we obtain 
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Result: u* 
Initialization: Set parameters as 0 = 7 = ar) and initial estimates as u9 = uz, 


au =0, vf of 0, v2, = 02, =0, vo = 02 =0, vf = of =0, ho =0, k° =0 and 
w° = 0. 
while convergence criterion not met, for k=0,1,... do 


Zk+1 Zk ~k _ a-ab 
1L-Akti = PYOX, (A ).|24.¢6) (h + o0(Fm(t& ))) PTOX,(A||.1|24.¢5) (2) a TeX 


2-BE +" = projj).|jcc<a (RY + o(Wa*)), (proj .<0(6)) GI) = 
min{max{s(i, j), —B}, }, Ving € {1,2,...,n}, 
3-akt1 = w* + o(OZOE + OF, v8, + Of0§ + Hf — Da*) 
d-ub+l = yk — (FA At 4 weet — De wet), 
ne 6 a 
Use formula Proxga,|.|,(@)(4, 9) = (1 ae eGietaot )a(i,j) for steps 5-8 
below: 
5-uy = Prox;|.|, (vf = 7(O,w**1)), 
6-uk +t = Prox;|.|, (vk, —7r(O4w*t1)), 
Tuk tt = Prox;|.|, (vk —7(O.w**?)), 
8-vitt = Prox,)), (vf! — 7r(O,w**")), 
9-ak+1 = Quk +h — ub, ort? = Qurtt — uf ot = Quhtt 


Syl yk att! = byt 


ko ogk+1 _ 
— UE,, Ve = 


vg 


k 
vy. 
end 


Algorithm 2: The algorithm for undersampled MRI reconstruction 


Uwe RM1xNo2 


kkt+1 kk Ww 0 0 0 0 vy E€ (R*)™2 x No 
hk+! | — prox, p= hk | +a Fm 0 0 0 0 vex € (Rt)§1*Ne2 
peti w* —D OF O% O§ O71 Ve € (R*)™2 x No 

VLE (R*)N1*Na2 


w + o(Oxvk + O*, vk, 


AES" = proxy aj i34.t6) (R* + o( Fin (u*))) Proxy |). 1134.t0)() = Tye 
KEY = projj).||.<o(k* + o(Wa")), (proj) <s(6)) (5) 


= min{max{s(i, 7), —S}, B}, Vi,7 € {Ui 2 ease 
wrt a< okt (Out + O%,vk, + Osvk + O*% vk — Du). 


4.1.2 The x-subproblems 


For x = (u, vp, Ves, Ve, v,)', according to Algorithm 1, we get, 
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yt uk w* Fe —p* 
of a 00" oO, | (RH 
yr tt | = prox,¢ OP |r| Oo Ab Des Akt 
gett uk 00 O wrt 
ue id 0 0 Of 
yk — 7(Utkett 4 rae _ D* w+) 
uf = T(O,w**") 
= prox.e wk, —7(O.,w*t1) , 
wk —7(O.0**") 
wk —7(O,wkt1) 


where 


D*u(n1,n2) = [ui(m1 — 1, n2) — ui (m1, n2)] + [ug(n1 — 1,2 — 1) — ue(n1, n2)] 
+[u3(n41 = 1, ng — 1) = u3(n1, N2)| 
+[ua(ny +1,n2—-1) — ua(n,n2)], (m1, n2) € A. 


Based on the above argument, we propose Algorithm 2, to solve undersampled 
MRI reconstruction problem. 


4.2 Image denoising formulation 


Let z € R%:*%2 be a gray scaled two-dimensional noisy image, and consider 
the following optimization problem: 


1 
min gllz — ull2 + ATVais(u). 


This problem is named denoising problem and can be rewritten with the 
following problem: 
: 1 2 * * * * 
min 5|[2— ul + logls + lotsa + loss + ey h 
s.t. Of (vt)(m1,n2) + OF (ve) (m1, M2) + OF (ve)(m1, n2) + OF (v+)(n1, n2) 


=Du(ni,n2), (m1,n2) € X. 


To solve this problem numerically, we need to determine functions F' and G 
as well as the linear operator K in (3). Assume x = (u, Ups Vass Ves vi)". Set 


iv 8 0 @ 1 
a (-n OF OF OF Oo) FY = 5|I¥: — 2ll2 + 5(0}(¥2 — Dz), 


Ge = |vpli + log|a + lodl + lop ls. 


Consequently 
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I —D* 
0 O; 
005),  Fty=|lmla+<zm >. 
00. 
0 O4 


Analogous to the arguments about Algorithm 2, we can obtain the similar 
procedure (Algorithm 3), to solve the image denoising problem. 


Result: u* 
Initialization: Set parameters as 0 = T= Fa and initial estimates as u° = z, 
~O 0, v9 ~0 0 0 


bi t Ct 0, ve, pan 0, ve Ue 0, vy Cie 0, n° 0 and we? = 0. 


while convergence criterion not met, for k=0,1,... do 

y-Aktl — PLOX,(1)).||24.tz) (h® + o(u*)), prox,(1)).1)24.¢2) (a) = a 
2-wktl = wk + o(OZvF + OF, HF, + OF0% + OF OF — Dak), 

B-yktl = uk = T(hktt ei D*w*t1), 


Use formula Proxga,|.|,(@)(4,9) = (1 mata laa) )a(i,j) for steps 4-7 


Ayktt = Prox,|.|, (vf - 7(O,w**1)), 
5-vtt = Prox), (vt, — T(O4w*t")), 
6-vet! = Prox,|.|, (vk — 7r(Oww*t1)), 
Tuy? = Prox,,), (vf — r(O,w**1)), 


t t ~ Up 
10-ok +1 = Quer — vk, 
1L-8+1 = Qyk+t _ yk, 
12-5477 = duttt — yk 

end 


Algorithm 3: The algorithm for denoising 


5 Numerical results 


5.1 Compressed MRI data reconstruction problem 


In this section, we run the proposed method on some test images and compare 
them with other well-known models. We choose a 256 x 256 in-vivo MR image, 
which is a T>-weighted axial brain scan. The image was chosen based on its 
variety of displayed feature and the fraction of FOV they occupy (about 46 
percent of the FOV). Based on these sizes, we chose to sample 10 percent 
of the k-space of the images. The random sampling masks were chosen as 
variable-density patterns, with more samples in the low-frequency area (in 
a neighborhood of the origin) and fewer ones towards the periphery of the 
k-space. We used the Sparse MRI MATLAB package [18] to generate the 
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masks. We compared the performance of our method with the baseline zero- 
filling solution, TV + wavelet model and the TGV model. The TV + wavelet 
model solves the unconstrained minimization problem 


. tl 
min 5||Fmu— O13 + ||Wulli + al|Vulli. 


This model was first introduced in [19] and then was further investigated in 
[20, 14, 16]. We followed the experiment setup proposed in these papers in 
our tests. In particular, we set a = 0.001 and 8 = 0.035. The TGV method 
has the form 


min 5 l|Fnu— O3 +ai||Vu — lle + allel. 
and was proposed in [15]. Again, in order to stay true to the original paper, 
we chose the parameters as in [15]. In particular, we set 4 = 8 x 10-°, 
a, = 1, and ap = 2. and the Chambolle—Pock algorithm is applied to 
them. In the proposed model (5), we also fix \ = 8 x 107° and 6 = 0.75 
across the simulations. The 4-tap Daubechies wavelet transform was used 
the experiment and was implemented through the Sparse MRI package. The 
algorithm ran for 500 iterations. 

We compare the quality of output images of each method based on SNR, 
SSIM, and HFEN indices. Let ug be the ground-truth image. The SNR is 


defined as 
SNR(u) := 20 log __ [tole 
||u — uol|e 


and the high SNR is related to low MSE. The SSIM index is defined as 


(2a + €1) (2025 19 C2) 
(ud + w2, + c1)(o% + oy2 + ¢2)’ 


SSIM(w) := 


where [Uy is the average of u, oy, is the standard deviation of u, uu, is the 
cross-covariance for u and ug, and c,; = (0.01 x 255)? and cy = (0.03 x 
255)? are regularizing constants. Another newly-proposed error metric is the 
HFEN (high frequency error norm), which aims at quantifying the quality of 
reconstructed edges and fine features [22]. This index is defined as 


||LoG(u) = LoG(uo)| lo 


HFEN = 
4) [[LoGuolle 


where LoG is a rotationally symmetric Laplacian of Gaussian filter with a 
kernel size of 15 x 15 pixels and standard deviation of 1.5 pixels. We note that 
the perfect reconstruction has SNR = oo, that SSIM = 1 and that HFEN 
= 0. We should also remark that none of these metrics perfectly quantify the 
visual perception of reconstructed images; although rarely the case, it might 


A four directions variational method for solving image processing problems 99 


0.3 i 
0.2 . 
0.1 - 
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Figure 1: T-weighted axial brain scan results. Top row: ground truth (left), sampling 
mask at rate %10 (middle), zero-filling solution (right). Middle row: reconstruction by 
TV+W (left), TGV (middle) and the proposed method (right). Bottom row: error maps 
for TV+W (left), TGV (middle) and the proposed method (right). 


just happen that between two reconstructions, one has a higher value in all 
of these metrics, yet the other is visually more faithful to the ground-truth. 


In Figure 1, the reconstruction results for the T>-weighted brain scan are 
presented. As expected, the zero-filling solution yields the worst reconstruc- 
tion with too many incoherent artifacts. The TV+W model removes many 
of these artifacts but still has low structural similarity to the ground truth, 
which is also easily observed in the reconstructed image. The TGV model 
significantly improves upon TV + wavelet model as can be perceived visually 
and quantitatively (Table 1). The proposed framework moderately improves 
upon the TGV solution as attested by all three error measures and the error 
maps (the third row in Figure 1), which shows fewer errors in [0, 0.1] interval. 
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Table 1: Quantitative comparison between various reconstructions for the T>-weighted 
brain scan. 


Method SNR [dB] SSIM HFEN 


Proposed 20.71 0.9046 0.1013 

TGV 20.07 0.9037 0.0995 

TV+W 17.58 0.5113 0.2259 

ZF 8.83 0.2733 0.7574 
22 
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Figure 2: SNR plots for the T2-weighted brain scan reconstructions. 


5.2 Denoising problem 


Consider denoising problem (4.2), assume that x is a clean grayscale image 
and z is a noisy image that is created by adding Gaussian noise with variance 
0.1 to the clean image x. In fact, we create an artificial noisy image by means 
of adding noise to the original image. Here we solve the denoising problem for 
the pirate test image and compare the results of the proposed model, with 
Condat [10], TGV [5], and Hosseini’s model [13]. Table 2 shows the denoised 
results. It can be seen that the proposed method has better performance 
in decreasing noise and preserving the details in comparison with the other 
methods. On the other hand, Figure 3 shows SSIM and PSNR values with 
respect to the steps. As a result, it can be seen that the proposed model has 
better values of SSIM and PSNR, at the solutions that are obtained in the 
convergence. 
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Table 2: Results of Pirate denoising problem. 


Reference 


Noisy 


Condat 


TGV 
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Figure 3: PSNR (left) and SSIM (right) values versus number of iterations for denoising 


problem. 


6 Conclusions 


In this paper, a new discretization of total variation functional was proposed. 
A four directions discrete gradient was used in the structure of the method. 
Moreover, some linear operators were inserted to the constraints. One of 
differences of the proposed method with the method in [13], is the new in- 
serted projecting operators to the corners of the pixels. The proposed method 
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was employed to formulate some image processing problems; reconstructing 
highly undersampled MR images, and image denoising. The Chambolle—Pock 
algorithm was applied for the numerical simulations. The wavelet transform 
was chosen in the formulation of compressed sensing MRI problem and it 
was compared with zero-filling solution, TGV, and TV+W. In addition, for 
the image denoising problem, the proposed model was compared with some 
state of the art methods. In both numerical experiments, the new methods 
illustrated better quantitative and qualitative results. 
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